%Script_for_determination_of_Velocity_Percentage_and_IndividualAnalysis(Similarity
%Coefficient?

%%VALID FOR RUNNING 5 EXPERIMENTS, FOR SPECIFIC TYPE OF REPEATED
%%EXPERIMENTS WHOSE STIMULI START AT SECOND 90
%%Yields Response Inde 1 (Max Stopping) Response Index 2 (Max Crawling
%%Backward) and change in response over the course of the xperiment fitted
%%to a linear plot (Slope)


%Step_1 Extract Experiment Data (MISC2B)%%%(OPTIONAL STEP)


%%Step 2 Input Parameters

prompt = 'How long is stimulus interval In Seconds? ';
    SoundLength = input(prompt);
prompt2 = 'How long is the experiment in seconds? ';
    ExperimentLength= input(prompt2);
prompt3 = 'How long is the total stimulus window (time from start of one stimulus to the start of the next) in seconds? ';
    WindowDuration= input(prompt3);
prompt4 = 'In which stimulus window does the stimulus start (1st=1, 2nd=2, 3rd=3, etc...?';
     ControlSequenceWindow=input(prompt4);
     NumExptper= ExperimentLength/WindowDuration;
     StimulusWindows= NumExptper-ControlSequenceWindow-1;

%%Step 3, Extract Behaviour Perccentages

   %%For this, it is important to identify what the behaviour we are
   %%looking for is exacttly. Based on observations and hypothetical
   %%conjecture, I would say the response of a larvae to a vibration is
   %%typically a stopping movment, a backwards crawl, and a turn, To identify correctlly what percentage of 
   %%larvae do which, it is important to recognize that each act happens in
   %%a logical order. i.e., a larvae cannot move backwards without first
   %%stopping, and a larvae cannot also turn without stopping. Thus, it
   %%will be useful to identify what percentage of larvae exhibit the
   %%following behaviours in the 10-second window of the stimulus:
   
       %%1: Stopping
       %%2: Stopping and Crawling Backwards
       %%3: Stopping and Crawling Backwards and Turning *** Hypothesis
       %%4: Stopping and Turning
       %%5: Crawling Normally
       
   %% On a final introductory note, it should be stated that all behaviour values should be based
    %% off of previously acquired control values for behaviour.
     SpeedSet= (1: ExperimentLength);
for i=1:length(eset.expt)  
    for j=1:length(eset.expt(i).track)  %%% 1 %%%
    times=0;
    HeadVec=0;
    HeadUnitVec=0;
    SpeedRun=0;
    SpeedRunVel=0;
        times = eset.expt(i).track(j).dq.eti;  %%% 2 %%%
        numPoints = length(times);
        % Positions:
        pos = eset.expt(i).track(j).getDerivedQuantity('sloc');   %%% 3 %%%
        xpos = pos(1,:);
        xpos = xpos*lengthPerPixel;
        ypos = pos(2,:);
        ypos = ypos*lengthPerPixel;
        % Combined Matrix:
        combined = [times;xpos;ypos];
        k=k+1;
        times= transpose (times);
        xpos= transpose (xpos);
        ypos= transpose (ypos);
        headpos=0;
        midpos=0;
        %%Extracting Unit Vector For Position of Head
        headpos=eset.expt(i).track(j).dq.shead;
        midpos=eset.expt(i).track(j).dq.smid;
        HeadVec=(headpos-midpos);
        for k=1:(size(HeadVec,2)-1)
            for w=1:2
                HeadUnitVec(w,k)= HeadVec(w,k)/(sqrt((HeadVec(1,k))^2+(HeadVec(2,k))^2));
            end
        end          
%%ExtractingSpeeds
for o=1:(size(times)-1)
    SpeedRun(o)= (sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))/(times(o+1)-times(o));
end
%%Idea for extracting velocity values by dotting the product of velocity
%%unit vector and position of head unit vector.
%%Final Speedrun is the linear one-dimensional velocity of the
%%forward movement of a larvae relative to the direction of the head
 VelocityVecx=0;
 VelocityVecy=0;
 VelocityVec=0;
 for o=1:(size(times)-1)     
    VelocityVecx(o) = [(xpos(o+1)-xpos(o))/(sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))];
    VelocityVecy(o) = [(ypos(o+1)-ypos(o))/(sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))];
   %% CosThetaFactor= dot(VelocityVec(l),HeadVec(l))
   %SpeedRun(o)= CosThetaFactor * (sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))/(times(o+1)-times(o));
 end
 %%VelocityVecx=transpose(VelocityVecx);
 %%VelocityVecy=transpose(VelocityVecy);
 VelocityVec=vertcat(VelocityVecx,VelocityVecy);
 %%Dotting Velocity and HeadOrientation
 VelocityRelativeToHead=dot(VelocityVec,HeadUnitVec);
 for o=1:(size(times)-1)    
     CosThetaFactor(o)=VelocityRelativeToHead(o)/ (sqrt(VelocityVec(1,o)^2+VelocityVec(2,o)^2)*sqrt(HeadUnitVec(1,o)^2+HeadUnitVec(2,o)^2)); 
     SpeedRunstop(o)= (sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))/(times(o+1)-times(o));
 end
%%Extracting Times Associated with Speeds
TimePoints=size(times);
times(TimePoints(1):end )=[];
%%Extracting Velocity SpeedRuns
for o=1:(size(times))
    SpeedRunVel(o) = SpeedRun(o)*CosThetaFactor(o);       
end
SpeedRun= transpose (SpeedRun);
SpeedRunVel= transpose (SpeedRunVel);
%%Interpolating Speeds and Times
InterpolatedSpeeds=0;
for SecondnumberStunted=1: ExperimentLength
    InterpolatedSpeeds(SecondnumberStunted)= interp1(times,SpeedRun,SecondnumberStunted);
end
InterpolatedSpeedsVel=0;

for SecondnumberStunted=1: ExperimentLength
    InterpolatedSpeedsVel(SecondnumberStunted)= interp1(times,SpeedRunVel,SecondnumberStunted);
end
%%RValueAnalysis
for b=0:NumExptper
    Mark1=0;
    MarkR1=0;
    for a=ControlSequenceWindow*WindowDuration:ExperimentLength
          if  b*WindowDuration<a &&(b*WindowDuration+SoundLength)>a && SpeedRunVel(a)< 1 
                Mark1=1;
          end
    end
    %%R1
    for a=ControlSequenceWindow*WindowDuration:ExperimentLength    
          if  (b+1)*WindowDuration<a &&((b+1)*WindowDuration+SoundLength)>a && SpeedRunVel(a)< 1 
                MarkR1=0.5;
          end
    end
          if (MarkR1+Mark1)==1.5
              R1Succeeds=R1Succeeds+1;
          end
          if (MarkR1+Mark1)==0.5
              R1Fails=R1Fails+1;
          end
          if (MarkR1+Mark1)==1
              R1Succeeds=R1Succeeds+1;  
          end
end     
Placeholder(j,:)=InterpolatedSpeeds;
PlaceholderVel(j,:)=InterpolatedSpeedsVel;
%%NEED TO FIX THIS, INEFFICIENT WAY OF EXTRACTING DATA. ESSENTIALLY, IN
%%ORDER TO EXTRACT THE DATA AS IT STANDS FROM MULTIPLE DIFFERENT
%%EXPERIMENTS, I NEED TO VERTCAT MULTIPLE ARRAYS INTO ONE LARGE ARRAY, Each
%%Extracted Array having its own set of larvae speeds and times. Each
%%Vertcated array represents a separate experiment set (i)
  %%The first Set is distinct from the rest
if i==1
    set1(j,:)=Placeholder(j,:);  
    set1Vel(j,:)=PlaceholderVel(j,:); 
end
if i==2
    set2(j,:)=Placeholder(j,:);
    set2Vel(j,:)=PlaceholderVel(j,:); 
end
if i==3
    set3(j,:)=Placeholder(j,:);
    set3Vel(j,:)=PlaceholderVel(j,:); 
end
if i==4
    set4(j,:)=Placeholder(j,:);
    set4Vel(j,:)=PlaceholderVel(j,:); 
end
    if i>4
    CombinedMatrix=vertcat(set1,set2,set3,set4,Placeholder);
    CombinedMatrixVel=vertcat(set1Vel,set2Vel,set3Vel,set4Vel,PlaceholderVel);
    end
end
end 
AveragedMatrix=mean(CombinedMatrix,'omitnan');
%%Finding Mean of first part
    for i=1:(85)
        ControlSpeednum(i)=AveragedMatrix(i);
    end
    %Extracting Speed from previous Set
    %%Define Stopping Threshold for Data Set based on Average Control Time 
    %5Extracted Previously
    %%Define Turn Theta Thresholds
       SpeedFactor = .1;
       StopThreshold = mean(ControlSpeednum)*SpeedFactor;
       TurnFactor1= (3.141/4);
       TurnFactor2= ((3*3.141)/4);
       TurnFactor3= ((5*3.141)/4);
       TurnFactor4= ((7*3.141)/4); 
       %%Defining Behavioural Parameters as null sets
       %%ERROR NEED TO FIX THIS; COULD NOT GET MAX LENGTH SIZE FOR SECOND
       %%ARGUMENT
      
    %%  LarvaeStop=NaN(length(eset.expt), 200 ,20);
      %%LarvaeStopBack=NaN(length(eset.expt), 200,20);
   %% LarvaeStopBackturn=NaN(length(eset.expt), 200,20);
   %% LarvaeStopturn=NaN(length(eset.expt), 200,20);
   %% LarvaeStopOP=NaN(length(eset.expt), 200,20);
for i=1:length(eset.expt)  
    for j=1:length(eset.expt(i).track)  %%% 1 %%%
        SpeedRunVel=0;
    times=0;
    HeadVec=0;
    HeadUnitVec=0;
    SpeedRun=0;
    SpeedRunVel=0;
        times = eset.expt(i).track(j).dq.eti;  %%% 2 %%%
        numPoints = length(times);
        % Positions:
        pos = eset.expt(i).track(j).getDerivedQuantity('sloc');   %%% 3 %%%
        xpos = pos(1,:);
        xpos = xpos*lengthPerPixel;
        ypos = pos(2,:);
        ypos = ypos*lengthPerPixel;
        % Combined Matrix:
        combined = [times;xpos;ypos];
        k=k+1;
        times= transpose (times);
        xpos= transpose (xpos);
        ypos= transpose (ypos);
        headpos=0;
        midpos=0;
        %%Extracting Unit Vector For Position of Head
        headpos=eset.expt(i).track(j).dq.shead;
        midpos=eset.expt(i).track(j).dq.smid;
        HeadVec=(headpos-midpos);
        for k=1:(size(HeadVec,2)-1)
            for w=1:2
                HeadUnitVec(w,k)= HeadVec(w,k)/(sqrt((HeadVec(1,k))^2+(HeadVec(2,k))^2));
            end
        end          
%%ExtractingSpeeds
for o=1:(size(times)-1)
    SpeedRun(o)= (sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))/(times(o+1)-times(o));
end
%%Idea for extracting velocity values by dotting the product of velocity
%%unit vector and position of head unit vector.
%%Final Speedrun is the linear one-dimensional velocity of the
%%forward movement of a larvae relative to the direction of the head
 VelocityVecx=0;
 VelocityVecy=0;
 VelocityVec=0;
 for o=1:(size(times)-1)     
    VelocityVecx(o) = [(xpos(o+1)-xpos(o))/(sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))];
    VelocityVecy(o) = [(ypos(o+1)-ypos(o))/(sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))];
 end
 %%VelocityVecx=transpose(VelocityVecx);
 %%VelocityVecy=transpose(VelocityVecy);
 VelocityVec=vertcat(VelocityVecx,VelocityVecy);
 %%Dotting Velocity and HeadOrientation
 VelocityRelativeToHead=dot(VelocityVec,HeadUnitVec);
 for o=1:(size(times)-1)    
     CosThetaFactor(o)=VelocityRelativeToHead(o)/ (sqrt(VelocityVec(1,o)^2+VelocityVec(2,o)^2)*sqrt(HeadUnitVec(1,o)^2+HeadUnitVec(2,o)^2)); 
     SpeedRunstop(o)= (sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))/(times(o+1)-times(o));
 end
%%Extracting Times Associated with Speeds
TimePoints=size(times);
times(TimePoints(1):end )=[];
%%Extracting Velocity SpeedRuns
SpeedRunStopCheck=NaN(size(times));
for o=1:(size(times))
    SpeedRunVel(o) = SpeedRun(o)*CosThetaFactor(o); 
    if  eset.expt(i).track(j).isrun(o) ==0
        SpeedRunStopCheck(o)=1; 
    else
        SpeedRunStopCheck(o)=0; 
    end
    
end
SpeedRun= transpose (SpeedRun);
SpeedRunVel= transpose (SpeedRunVel);
SpeedRunStopCheck=transpose(SpeedRunStopCheck);
for SecondnumberStunted=1: ExperimentLength
InterpolatedSpeedsStopCheck(SecondnumberStunted)= interp1(times,SpeedRunStopCheck,SecondnumberStunted);
end
InterpolatedSpeedsVel=0;
for SecondnumberStunted=1: ExperimentLength
InterpolatedSpeedsVel(SecondnumberStunted)= interp1(times,SpeedRunVel,SecondnumberStunted);
if InterpolatedSpeedsVel(SecondnumberStunted)<0
    Checker(SecondnumberStunted)=1;
else
    Checker(SecondnumberStunted)=0;
end

end
PlaceholderStopCheck(j,:)=InterpolatedSpeedsStopCheck;
PlaceholderVel(j,:)=(Checker);

%%NEED TO FIX THIS, INEFFICIENT WAY OF EXTRACTING DATA. ESSENTIALLY, IN
%%ORDER TO EXTRACT THE DATA AS IT STANDS FROM MULTIPLE DIFFERENT
%%EXPERIMENTS, I NEED TO VERTCAT MULTIPLE ARRAYS INTO ONE LARGE ARRAY, Each
%%Extracted Array having its own set of larvae speeds and times. Each
%%Vertcated array represents a separate experiment set (i)
  %%The first Set is distinct from the rest
if i==1
    set1(j,:)=PlaceholderStopCheck(j,:);  
    set1Vel(j,:)=PlaceholderVel(j,:); 
end
if i==2
    set2(j,:)=PlaceholderStopCheck(j,:);
    set2Vel(j,:)=PlaceholderVel(j,:); 
end
if i==3
   set3(j,:)=PlaceholderStopCheck(j,:);
    set3Vel(j,:)=PlaceholderVel(j,:); 
end
if i==4
    set4(j,:)=PlaceholderStopCheck(j,:);
    set4Vel(j,:)=PlaceholderVel(j,:); 
end
    if i>4
    CombinedMatrixStopCheck=vertcat(set1,set2,set3,set4,PlaceholderStopCheck);
    CombinedMatrixVel=vertcat(set1Vel,set2Vel,set3Vel,set4Vel,PlaceholderVel);
    end
%%%Extracting Stop Behaviour (u represents the numbered 30-sec
%%window)

%for o=1:(size(times)-1)
    %%Filling in non-stop holes in NAN array first
    %for u=1:20
    %if ((times(o) < (30*u)+10))  && (times(o) > 30*u)
        %LarvaeStop(i,j,u)=0;
         %LarvaeStopBack(i,j,u)=0;
         % LarvaeStopBackturn(i,j,u)=0;
          %LarvaeStopturn(i,j,u)=0;
          %LarvaeStopOP(i,j,u)=0;
    %else     
    %end
    %end
    %%Extracting Stops
%for u=1:20
   % if ( eset.expt(i).track(j).isrun(o) ==0 && (times(o) < (30*u)+10))  && (times(o) > 30*u)
     % LarvaeStop(i,j,u)=1;
       %LarvaeStopOP(i,j,u)= LarvaeStopOP(i,j,u)+1;
   % else   
      %  mello=2;
    %end
 %end
%%Extracting Stop and Back Behaviour 
%for u=1:20
    %if (eset.expt(i).track(j).isrun(o) ==0 && (times(o) < (30*u)+10))  && (times(o) > 30*u && SpeedRunVel(o)< 0 )
      % LarvaeStopBack(i,j,u)=1;
    %else      
   % end
%end
%%%Extracting and Stop and Back and Turn
%for u=1:20
   %if (eset.expt(i).track(j).isrun(o) ==0 && (times(o) < (30*u)+10))  && (times(o) > 30*u && SpeedRunVel(o)< 0 && abs(acos(CosThetaFactor(o))) > TurnFactor1 && abs(acos(CosThetaFactor(o))) < TurnFactor2 )
     %  LarvaeStopBackturn(i,j,u)=1;
    %else 
       % end 
%%Extracting and Stop  and Turn only
%for u=1:20
   % if (eset.expt(i).track(j).isrun(o) ==0 && (times(o) < (30*u)+10))  && (times(o) > 30*u && abs(acos(CosThetaFactor(o))) > TurnFactor1 && abs(acos(CosThetaFactor(o))) < TurnFactor2)
      %  LarvaeStopturn(i,j,u)=1;
   % else 
      %  end
    %end
    
end
end
%end
%end 
%%Final Behavioural Compilation
%LarvaeStopX=mean(LarvaeStop,2,'omitnan');
%StopPercentageY=mean(LarvaeStopX,1,'omitnan');
%StopPercentage=mean(StopPercentageY,3,'omitnan')
%LarvaeStopOPX=mean(LarvaeStopOP,2,'omitnan');
%StopPercentageOPY=mean(LarvaeStopOPX,1,'omitnan');
%StopPercentageOP=mean(StopPercentageOPY,3,'omitnan')
%OpportunityPercentage=(StopPercentageOP/10)
%ContinuePercentage= 1-StopPercentage
%LarvaeStopBackX=mean(LarvaeStopBack,2,'omitnan');
%StopBackPercentageY=mean(LarvaeStopBackX,1,'omitnan');
%StopBackPercentage=mean(StopBackPercentageY,3,'omitnan')
%LarvaeStopBackturnX=mean(LarvaeStopBackturn,2,'omitnan');
%StopBackturnPercentageY=mean(LarvaeStopBackturnX,1);
%StopBackturnPercentage=mean(StopBackturnPercentageY,3,'omitnan')
%LarvaeStopturnX=mean(LarvaeStopturn,2,'omitnan');
%StopturnPercentageY=mean(LarvaeStopturnX,1,'omitnan');
%StopturnPercentage=mean(StopturnPercentageY,3,'omitnan')
   % end 
%end
%%Plotting overalL Velocity throughout 600 Seconds
%%Extracting Average Speedvalues from combined matrix for each time point
%%which does not exhibit a NAN,taking the mean of the Column row vectors at
%%each point
AveragedMatrixStopCheck=mean(CombinedMatrixStopCheck,'omitnan');
AveragedMatrixVel=((mean(CombinedMatrixVel,'omitnan')));
subplot(4,1,1)
hold on
plot(SpeedSet,AveragedMatrixVel,'b')
title('Proportion of Larvae Crawling Backward')
axis([0 ExperimentLength min(AveragedMatrixVel) max(AveragedMatrixVel)])
xlabel('Seconds Passed')
ylabel('Proportion of Larvae Crawling Backward')
hold off

 ControlSpeed2=mean(ControlSpeednum);

 subplot(4,1,2)
 hold on
plot(SpeedSet,AveragedMatrixStopCheck,'b')
title('Proportion of Larvae Stopped')
axis([0 ExperimentLength 0 1])
xlabel('Seconds Passed')
ylabel('Proportion of Larvae Stopped')
hold off

    %%Still need to: Extract Standard Errors    
%%Step 4 Integral Analysis

%%Plotting Averaged Velocity over Condensed  window starting from
%%Stimulus Windows
Window=zeros(1,WindowDuration);
for SecondnumberStunted=1:WindowDuration
            for p=1:(NumExptper+1-ControlSequenceWindow)
                for j=1:ExperimentLength
                    if ((WindowDuration*(p-1)+SecondnumberStunted-WindowDuration+90)==j)
                         Window(SecondnumberStunted)=Window(SecondnumberStunted)+AveragedMatrixVel(j);
                         
                    end
                end
            end
end   
 Window= Window/(StimulusWindows);
 subplot(4,1,3)
 hold on
 fill(WindowDuration,max(Window),'r')
colormap(cool)
bar(SoundLength/2,max(Window),SoundLength) 
 TimeofWIndow= (1:WindowDuration);
 plot(TimeofWIndow,Window)
  hold off
   subplot(4,1,4)
   hold on
  Maxnum=zeros(1,NumExptper);
  IntVel=zeros(1,NumExptper);
for TimeWindows=1:NumExptper-1
    for i=0:WindowDuration-1
        for j=1:ExperimentLength-4
        if ((TimeWindows)*WindowDuration+i+90)==j && AveragedMatrixVel(j) > Maxnum(TimeWindows)
            Maxnum(TimeWindows)= AveragedMatrixVel(j);
            IntVel(TimeWindows)= (AveragedMatrixVel(j)+AveragedMatrixVel(j+1)+AveragedMatrixVel(j+2))/(3);
        end
        end
    end
end
ix =(Maxnum==0);
Maxnum(ix)=[];
iy =(IntVel==0);
IntVel(iy)=[];
Experimentstotal=(1:length(IntVel));
Experimentstotal2=(1:length(Maxnum));
Maxnum2=zeros(1,NumExptper);
IntVel2=zeros(1,NumExptper);
Q1=polyfit(Experimentstotal,IntVel,1);
for TimeWindows=1:NumExptper-1
    for i=10:WindowDuration-10
        for j=1:ExperimentLength
        if ((TimeWindows)*WindowDuration+i+90)==j && AveragedMatrixVel(j) > Maxnum2(TimeWindows)
            Maxnum2(TimeWindows)= AveragedMatrixVel(j);
            IntVel2(TimeWindows)= (AveragedMatrixVel(j)+AveragedMatrixVel(j+2)+AveragedMatrixVel(j+1))/(3);
            
        end
        end
    end
end
iq =(Maxnum2==0);
Maxnum2(iq)=[];
it =(IntVel2==0);
IntVel2(it)=[];
Experimentstotals=(1:length(IntVel2));
Experimentstotal2s=(1:length(Maxnum2));
Q2=polyfit(Experimentstotals,IntVel2,1);
x1=linspace(0,StimulusWindows,100);
Fitted1=polyval(Q1,x1);
Fitted2=polyval(Q2,x1);
plot(Experimentstotal,IntVel,'b',Experimentstotal2,Maxnum,'r',Experimentstotals,IntVel2,'g',Experimentstotal2s,Maxnum2,'m')
plot(x1,Fitted1)
plot(x1,Fitted2)
hold off
Slope=Q1(1)
ResponseIndex2=max(AveragedMatrixVel)
ResponseIndex1=max(AveragedMatrixStopCheck)
R1=(R1Fails)/(R1Succeeds)
path of the larvae and show behaviours at different
       %%Time points  (.mov files might be useful for this analysis)
    %%Show speed Graphs for individual Larvae
   
%%ExtractingSpeeds
%%Finding whether a larvae stopped or didn't stop per stimulus number

